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Abstract 

Battery  management  systems  in  hybrid  electric  vehicle  battery  packs  must  estimate  values  descriptive  of  the  pack’s  present  operating 
condition.  These  include:  battery  state  of  charge,  power  fade,  capacity  fade,  and  instantaneous  available  power.  The  estimation  mechanism 
must  adapt  to  changing  cell  characteristics  as  cells  age  and  therefore  provide  accurate  estimates  over  the  lifetime  of  the  pack. 

In  a  series  of  three  papers,  we  propose  a  method,  based  on  extended  Kalman  filtering  (EKF),  that  is  able  to  accomplish  these  goals  on 
a  lithium  ion  polymer  battery  pack.  We  expect  that  it  will  also  work  well  on  other  battery  chemistries.  These  papers  cover  the  required 
mathematical  background,  cell  modeling  and  system  identification  requirements,  and  the  final  solution,  together  with  results. 

In  order  to  use  EKF  to  estimate  the  desired  quantities,  we  first  require  a  mathematical  model  that  can  accurately  capture  the  dynamics  of 
a  cell.  In  this  paper  we  “evolve”  a  suitable  model  from  one  that  is  very  primitive  to  one  that  is  more  advanced  and  works  well  in  practice. 
The  final  model  includes  terms  that  describe  the  dynamic  contributions  due  to  open-circuit  voltage,  ohmic  loss,  polarization  time  constants, 
electro-chemical  hysteresis,  and  the  effects  of  temperature.  We  also  give  a  means,  based  on  EKF,  whereby  the  constant  model  parameters 
may  be  determined  from  cell  test  data.  Results  are  presented  that  demonstrate  it  is  possible  to  achieve  root-mean- squared  modeling  error 
smaller  than  the  level  of  quantization  error  expected  in  an  implementation. 
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1.  Introduction 

This  paper  is  the  second  in  a  series  of  three  that  de¬ 
scribe  advanced  algorithms  for  a  battery  management  sys¬ 
tem  (BMS)  for  hybrid  electric  vehicle  (HEV)  application. 
This  BMS  is  able  to  estimate  battery  state  of  charge  (SOC), 
power  fade,  capacity  fade  and  instantaneous  available  power, 
and  is  able  to  adapt  to  changing  cell  characteristics  over  time 
as  the  cells  in  the  battery  pack  age.  The  algorithms  have 
been  implemented  on  a  lithium-ion  polymer  battery  (LiPB) 
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pack,  but  we  expect  them  to  work  well  for  other  battery 
chemistries. 

The  method  that  we  use  to  estimate  these  parameters  is 
based  on  Kalman  filter  theory.  (There  have  been  other  re¬ 
ported  methods  for  SOC  estimation  that  use  Kalman  filter¬ 
ing  [1,2],  but  the  method  in  this  series  of  papers  expands  on 
these  results  and  also  differs  in  some  important  respects,  as 
will  be  outlined  later.)  Kalman  filters  are  an  intelligent — and 
sometimes  optimal — means  for  estimating  the  state  of  a  dy¬ 
namic  system.  By  modeling  our  battery  system  to  include 
the  wanted  unknown  quantities  in  the  “state”,  we  may  use 
the  Kalman  filter  to  estimate  their  values.  An  additional  ben¬ 
efit  of  the  Kalman  filter  is  that  it  automatically  provides 
dynamic  error-bounds  on  these  estimates  as  well.  We  ex¬ 
ploit  this  fact  to  give  aggressive  performance  from  our  bat¬ 
tery  pack,  without  fear  of  causing  damage  by  overcharge  or 
overdischarge. 
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The  first  paper  [3]  is  an  introduction  to  the  problem.  It 
describes  the  HEV  environment  and  the  requirement  spec¬ 
ifications  for  a  BMS.  The  remainder  of  the  paper  is  a  brief 
tutorial  on  the  Kalman  filter  theory  necessary  to  grasp  the 
content  of  the  remaining  papers;  additionally,  a  nonlinear 
extension  called  the  “extended  Kalman  filter”  (EKF)  is 
discussed. 

This  second  paper  describes  some  mathematical  cell  mod¬ 
els  that  may  be  used  with  this  method.  The  HEV  application 
is  a  very  harsh  environment,  with  rate  requirements  up  to  and 
exceeding  ±20C  and  very  dynamic  rate  profiles.  This  is  in 
contrast  to  relatively  benign  portable-electronic  applications 
with  constant  power  output  and  fractional  C  rates.  Methods 
for  estimating  SOC  that  work  well  in  portable-electronic  de¬ 
vices  may  not  work  well  in  the  HEV  application.  If  precise 
SOC  estimation  is  required  by  the  HEV,  then  a  very  accurate 
cell  model  is  necessary. 

Results  of  lab  tests  on  physical  cells  are  presented  and 
compared  with  model  prediction.  The  best  modeling  results 
obtained  to  date  are  so  precise  that  the  root-mean- squared 
(RMS)  estimation  error  is  less  than  the  quantization  noise 
floor  expected  in  our  battery  management  system  design. 
More  importantly,  the  model  allows  very  precise  SOC 
estimation,  therefore  allowing  the  vehicle  controller  to  con¬ 
fidently  use  the  battery  pack’s  full  operating  range  without 
fear  of  over-  or  under-charging  cells.  This  paper  also  gives 
an  overview  of  other  modeling  methods  in  the  literature 
and  shows  how  an  EKF  may  be  used  to  adaptively  identify 
unknown  parameters  in  a  cell  model,  in  real  time,  given  test 
data. 

The  third  paper  [4]  covers  the  real-time  parameter  estima¬ 
tion  problem;  namely,  how  to  dynamically  estimate  SOC, 
power  fade,  capacity  fade,  available  power  and  so  forth.  An 
EKF  is  used  in  conjunction  with  the  cell  model.  The  cell 
model  may  be  fixed,  or  may  itself  have  adaptable  parame¬ 
ters  so  that  the  model  tracks  cell  aging  effects.  Details  for  a 
practical  implementation  are  discussed. 

We  now  proceed  by  briefly  reviewing  cell  models  in  the 
literature  that  have  been  proposed  for  SOC  estimation.  We 
explain  why  these  do  not  meet  the  requirements  presented 
in  [3].  Several  models  from  Refs.  [5,6]  do  meet  the  require¬ 
ments,  and  they  are  described  in  detail  here,  together  with 
some  new  models  and  results.  A  method  for  identifying 
model  parameters  using  an  extended  Kalman  filter  is  pre¬ 
sented,  followed  by  conclusions. 

2.  Standard  cell-modeling  methods  for  SOC  estimation 

The  literature  documents  a  number  of  cell-modeling  meth¬ 
ods  for  SOC  estimation.  An  excellent  summary,  in  greater 
detail  than  can  be  presented  here,  may  be  found  in  refer¬ 
ence  [7].  Here,  we  investigate  to  see  whether  any  of  these 
methods  meets  our  needs.  Recall  that  our  application  is  to 
model  cell  dynamics  for  the  purpose  of  SOC  estimation  in 
an  HEV  battery  pack. 


For  this  application,  the  cell  model  must  be  accurate  for  all 
operating  conditions.  These  include:  very  high  rates  (up  to 
about  ±20C,  far  exceeding  the  low  rates  considered  by  many 
papers  in  the  literature  for  portable  electronic  applications), 
temperature  variation  in  the  automotive  range  of  —30  to 
50  °C,  very  dynamic  rates  (unlike  the  more  benign  portable 
electronic  and  battery  electric  vehicle  application).  Charging 
must  be  accounted  for  in  the  model. 

We  also  require  non-invasive  methods  using  only  readily 
available  signals.  This  requirement  is  imposed  by  the  HEV 
environment  where  the  BMS  has  no  direct  control  over  cur¬ 
rent  and  voltage  experienced  by  the  battery  pack — these  are 
in  the  domain  of  the  vehicle  controller  and  inverter.  We  must 
rely  on  such  measurements  as  instantaneous  cell  terminal 
voltage,  cell  current  and  cell  external  temperature. 

Our  cell  chemistry  also  limits  the  range  of  approaches  we 
might  consider.  Techniques  specific  to  lead-acid  chemistries, 
for  example,  are  not  appropriate  for  LiPB  cells. 

2.7.  Laboratory  and  chemistry -dependent  methods 

Several  methods  for  direct  SOC  estimation  simply  cannot 
be  used  in  our  application: 

1.  A  laboratory  method  for  determining  SOC  is  to 
completely  discharge  a  cell,  recording  discharged 
ampere-hours,  to  determine  its  present  remaining  capac¬ 
ity.  This  is  the  most  accurate  SOC  measurement  tech¬ 
nique,  but  is  impractical  in  HEV  as  the  battery  energy 
is  wasted  by  the  test,  and  the  test  cannot  dynamically 
estimate  SOC. 

2.  Chemistry-dependent  methods  for  other  chemistries, 
such  as  Coup  de  Fouet  measurement,  or  measurement 
of  electrolyte  physical  properties  for  lead-acid  batteries, 
are  all  inappropriate  (as  our  application  uses  LiPB  cells). 

3.  Open-circuit  voltage  (OCV)  measurements:  If  the  cell  is 
allowed  to  rest  for  a  long  period,  its  terminal  voltage  de¬ 
cays  to  OCV,  and  OCV  may  be  used  to  infer  SOC  (via 
lookup  table,  for  example).  However,  long  periods  (some¬ 
times  hours)  of  battery  inactivity  must  occur  before  the 
terminal  voltage  approaches  OCV.  This  method  may  not 
be  used  for  dynamic  SOC  estimation.  (Other  complica¬ 
tions  with  this  method  include  the  dependence  of  OCV 
on  temperature,  and  presence  of  terminal  voltage  hys¬ 
teresis,  especially  at  low  temperatures.) 

2.2.  Electro -chemical  modeling 

One  approach  to  modeling  cell  electrical  dynamics  is 
to  carefully  consider,  at  the  molecular  level,  the  various 
processes  that  occur  within  the  cell.  Accurate  terminal 
voltage  prediction  may  be  achieved  by  these  models  (see 
Ref.  [8],  for  example).  However,  it  would  be  difficult  (if 
possible)  to  measure  the  many  required  physical  param¬ 
eters  on  a  cell-by-cell  basis  in  a  high-volume  consumer 
product.  We  have  not  pursued  this  approach,  although 
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we  employ  many  of  the  macroscopic  concepts  in  our 
models. 

2.3.  Impedance  spectroscopy 

Another  broad  category  of  cell  modeling  involves  mea¬ 
suring  cell  impedances  over  a  wide  range  of  ac  frequencies 
at  different  states  of  charge  [9-13].  Values  of  the  model 
parameters  are  found  by  least-squares  fitting  to  measured 
impedance  values.  SOC  may  be  indirectly  inferred  by  mea¬ 
suring  present  cell  impedance  and  correlating  them  with 
known  impedances  at  various  SOC  levels.  We  must  also  dis¬ 
count  this  method  for  our  application,  as  we  have  no  direct 
method  to  inject  signals  into  cells  to  measure  impedances. 
We  rely  on  the  vehicle  to  generate  and  dissipate  all  energy 
flowing  through  the  battery  pack. 

2.4.  Circuit  models 

A  number  of  papers  present  equivalent  circuit  models  of 
cells  [14-17].  Typically,  a  high- valued  capacitor  or  voltage 
source  is  used  to  represent  the  open-circuit  voltage.  The  re¬ 
mainder  of  the  circuit  models  the  cell’s  internal  resistance 
and  more  dynamic  effects  such  as  terminal  voltage  relax¬ 
ation.  From  the  OCV  estimate,  SOC  may  be  inferred  via 
table  lookup.  Both  linear-  and  nonlinear-circuit  models  may 
be  used.  Our  model  has  many  similarities  to  a  circuit  model, 
except  that  our  fundamental  “state”  is  SOC,  not  OCV. 

2.5.  Coulomb  counting 

The  final  method  discussed  in  the  literature  involves  SOC 
estimation  directly  via  Coulomb  counting.  This  may  be  done 
“open-loop”,  which  is  very  sensitive  to  current  measurement 
error,  or  “closed-loop”  which  can  be  much  more  accurate. 
The  feedback  mechanism  may  be  empirically  designed  [18] 
or  may  use  a  mathematically  optimized  approach  such  as 
the  Kalman  filtering  method  [1,2]  to  generate  the  feedback. 
All  Kalman  filtering-based  methods  by  other  authors  in  the 
literature  (with  which  we  are  familiar)  use  a  circuit  model  of 
the  cell  with  voltage  sources  and  capacitor  voltages  repre¬ 
senting  OCV  and  relaxation  effects.  OCV  may  be  estimated 
and  SOC  inferred  from  OCV. 

Our  approach  is  also  based  on  the  Kalman  filtering 
method,  but  the  fundamental  aspect  of  our  model  that  sets  it 
apart  from  those  reported  in  the  literature  is  that  SOC  itself 
is  required  to  be  a  state  of  the  system.  The  direct  benefit  of 
this  approach  is  that  the  Kalman  filter  automatically  gives 
a  dynamic  estimate  of  the  SOC  and  its  uncertainty  (this 
is  discussed  in  greater  detail  in  Ref.  [4]).  That  is,  instead 
of  reporting  the  SOC  to  the  vehicle  controller  (at  some 
point  in  time)  to  be  “about  55%”,  the  algorithm  is  able  to 
report  that  the  SOC  is  55  zb  3%,  for  example.  This  allows 
the  vehicle  controller  to  confidently  use  the  battery  pack’s 
full  operating  range  without  fear  of  over-  or  under-charging 
cells. 


3.  An  evolution  of  cell  model  structures 

In  order  to  use  Kalman-based  methods  for  a  battery 
management  system,  we  must  first  have  a  cell  model  in  a 
discrete-time  state-space  form.  Specifically,  we  assume  the 
form 


Xk+ 1  =  f(Xk,  Uk)  +  Wk, 

(1) 

yk  =  g(xk,  uk)  +  vk, 

(2) 

where  Xk  is  the  system  state  vector  at  discrete-time  index  k , 
where  the  “state”  of  a  system  comprises  in  summary  form 
the  total  effect  of  past  inputs  on  the  system  operation  so 
that  the  present  output  may  be  predicted  solely  as  a  function 
of  the  state  and  present  input.  Values  of  past  inputs  are  not 
required.  The  vector  Uk  is  the  measured  exogenous  system 
input  at  time  k  and  Wk  is  unmeasured  “process  noise”  that 
affects  the  system  state.  The  system  output  is  y k  and  is 
measurement  noise,  which  does  not  affect  the  system  state. 
Equation  (1)  is  called  the  “state  equation”,  (2)  is  called  the 
“output  equation”,  and  /(•,  •)  and  g(-,  •)  are  (possibly  non¬ 
linear)  functions,  specified  by  the  particular  cell  model  used. 
All  of  the  system  dynamics  are  represented  in  (1).  Equation 
(2)  is  a  static  relationship.  In  the  models  to  follow,  Wk  is 
used  to  account  for  current- sensor  error  and  inaccuracy  of 
the  state  equation,  and  is  used  to  account  for  voltage  sen¬ 
sor  error  and  inaccuracy  of  the  output  equation. 

In  the  case  where  we  wish  to  model  a  cell’s  dynamics 
using  (1)  and  (2),  the  vector  Uk  contains  the  instantaneous 
cell  current  4-  It  may  also  contain  the  cell  temperature  7^, 
an  estimate  of  the  cell’s  capacity  C,  and/or  an  estimate  of 
the  cell’s  internal  resistance  Rk ,  for  example.  The  system 
output  is  typically  a  scalar  but  may  be  vector  valued  as  well. 
Here  we  consider  the  output  to  be  the  cell’s  loaded  terminal 
voltage  (not  its  at-rest  OCV).  Our  method  constrains  the 
state  vector  Xk  to  include  SOC  as  one  component. 

There  are  many  possible  candidates  for  (1)  and  (2),  and  for 
the  choice  of  Xk  and  u^.  Here,  we  describe  the  development 
of  some  modeling  equations  in  an  evolutionary  sense.  That 
is,  we  start  with  a  very  simple  model,  and  gradually  add 
complexity  to  better  represent  the  true  cell  dynamics.  In 
order  to  justify  the  changes  in  the  model,  we  compare  model 
dynamics  to  cell  dynamics  based  on  data  collected  from  cell 
tests.  We  first  describe  the  cell  tests,  and  then  develop  the 
model  structures. 

3.1.  Cell  tests  for  model  fitting 

In  order  to  compare  the  abilities  of  the  proposed  models  to 
capture  a  cell’s  dynamics,  we  gathered  data  from  a  prototype 
LiPB  cell.  The  cell  comprises  a  LfiVh^CU  cathode,  an  artifi¬ 
cial  graphite  anode,  is  designed  for  high-power  applications, 
has  a  nominal  capacity  of  7.5  Ah  and  a  nominal  voltage  of 
3.8  V.  For  the  tests,  we  used  a  Tenney  thermal  chamber  set  at 
25 °C  and  an  Arbin  BT2000  cell  cycler.  Each  channel  of  the 
Arbin  was  capable  of  20  A  current,  and  ten  channels  were 
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SOC  and  current  as  a  function  of  time  during  discharge 


SOC  and  current  as  a  function  of  time  during  charge 


Fig.  1.  Plots  showing  SOC  vs.  time  and  rate  vs.  time  for  pulsed-current  cell  tests.  Discharge  portion  of  test  is  shown  in  (a);  charge  portion  of  test  is 
shown  in  (b).  Dark  line  is  SOC,  gray  line  is  current. 


connected  in  parallel  to  achieve  currents  of  up  to  200  A.  The 
cycler’s  voltage  measurement  accuracy  was  ±5  mV  and  its 
current  measurement  accuracy  was  ±200  mA. 

Two  types  of  cell  tests  were  performed  for  the  work  re¬ 
ported  in  this  paper.  The  first  type  comprised  a  sequence  of 
constant-current  discharge  pulses  and  rests  followed  by  a  se¬ 
quence  of  constant-current  charge  pulses  and  rests.  The  cell 
started  fully  charged  before  the  test  began.  Discharge  current 
pulses  from  150  down  to  1  A,  and  charge  pulses  from  150 
down  to  1  A  were  used.  The  current  and  SOC  profiles  for 
this  test  are  shown  in  Fig.  1(a)  and  (b).  Frame  (a)  shows  the 
discharge  portion  of  the  test  and  frame  (b)  shows  the  charge 
portion  of  the  test.  Data  points  (including  voltage,  current, 
ampere-hours  discharged  and  ampere-hours  charged)  were 
collected  once  per  second. 

The  second  test  was  a  sequence  of  16  urban  dynamome¬ 
ter  driving  schedule  (UDDS)  cycles,  separated  by  40  A  dis¬ 
charge  pulses  and  5  min  rests,  and  spread  over  the  90-10% 
SOC  range.  The  SOC  as  a  function  of  time  is  plotted  in 
Fig.  2(a),  and  rate  as  a  function  of  time  for  one  of  the  UDDS 
cycles  is  plotted  in  Fig.  2(b).  We  see  that  SOC  increases  by 
about  5%  during  each  UDDS  cycle,  but  is  brought  down 
about  10%  during  each  discharge  between  cycles.  The  en¬ 
tire  operating  range  for  these  cells  (10-90%  SOC)  is  excited 
during  the  cell  test. 


SOC  as  a  function  of  time 


The  data  was  used  to  identify  parameters  of  the  cell  mod¬ 
els  to  be  described  in  the  next  sections.  The  goal  is  to  have 
the  cell  model  output  resemble  the  cell  terminal  voltage  un¬ 
der  load  as  closely  as  possible,  at  all  times,  when  the  cell 
model  input  is  equal  to  the  cell  current.  Model  fit  was  judged 
by  comparing  root-mean-squared  estimation  error  (estima¬ 
tion  error  equals  cell  voltage  minus  model  voltage)  over  the 
portions  of  the  cell  tests  where  SOC  was  between  5  and 
95%.  Model  error  outside  that  SOC  range  was  not  consid¬ 
ered  as  the  HEV  pack  operation  design  limits  are  10-90% 
SOC. 

3.2.  SOC  as  a  state-vector  component 

The  one  constant  requirement  for  the  cell  model  is  that 
we  constrain  SOC,  denoted  as  Zk,  to  be  a  member  of  the 
state  vector  Xk.  To  be  careful,  we  give  a  list  of  definitions 
culminating  in  our  understood  definition  of  SOC. 

•  Definition'.  A  cell  is  fully  charged  when  its  voltage  reaches 
v  —  v\i  after  being  charged  at  infinitesimal  current  levels. 
Here,  we  use  v h  =  4.2  V  at  room  temperature  (25  °C). 

•  Definition'.  A  cell  is  fully  discharged  when  its  voltage 
reaches  v  =  v\  after  being  drained  at  infinitesimal  current 
levels.  Here,  we  use  v\  =  3.0  V  at  room  temperature. 


Current  for  one  UDDS  cycle  (zoom) 


Fig.  2.  Plots  showing  SOC  vs.  time  and  rate  vs.  time  for  UDDS  cell  tests.  SOC  is  shown  in  (a);  rate  for  one  UDDS  cycle  is  shown  in  (b). 
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•  Definition:  The  capacity  of  a  cell  C  is  the  maximum  num¬ 
ber  of  ampere-hours  that  can  be  drawn  from  the  cell  be¬ 
fore  it  is  fully  discharged,  at  room  temperature,  starting 
with  the  cell  fully  charged. 

•  Definition :  The  nominal  capacity  Cn  of  the  cell  is  the 
number  of  ampere-hours  that  can  be  drawn  from  the  cell 
at  room  temperature  at  the  C/30  rate,  starting  with  the  cell 
fully  charged. 

•  Definition :  The  SOC  of  the  cell  is  the  ratio  of  the  remain¬ 
ing  capacity  to  the  nominal  capacity  of  the  cell,  where  the 
remaining  capacity  is  the  number  of  ampere-hours  that 
can  be  drawn  from  the  cell  at  room  temperature  at  the 
C/30  rate  before  it  is  fully  discharged. 

With  these  definitions  in  place,  we  can  then  investigate 
some  mathematical  relations  involving  SOC.  Particularly, 

C  rjfii  r) 

z(t)  =  z( 0)  -  /  -!l±±  dr,  (3) 

JO  Cn 

where  z(t)  is  the  cell  SOC,  i(t)  is  instantaneous  cell  current 
(assumed  positive  for  discharge,  negative  for  charge),  and 
Cn  is  the  cell  nominal  capacity.  Cell  Coulombic  efficiency 
rji  is  rji  =  1  for  discharge,  and  rji  =  0  <  1  for  charge. 

Using  a  rectangular  approximation  for  integration  and  a 
“suitably  small”  sampling  period  At,  a  discrete-time  approx¬ 
imate  recurrence  may  then  be  written  as 


In  these  models,  yk  is  the  cell  terminal  voltage,  R  is  the 
cell  internal  resistance  (different  values  may  be  used  for 
charge/discharge  and  at  different  SOC  levels  if  desired),  Kt  is 
the  polarization  resistance  and  K\,  and  K 3  are  constants 
chosen  to  make  the  model  fit  the  data  well.  All  of  the  terms  of 
these  models  may  be  collected  to  make  a  “combined  model” 
that  performs  better  than  any  of  the  individual  models  alone. 
This  model  is 


Zk+ 1  —  Zk 


rfAfi 

Cn 


Ik, 


K 1 

yk  =  K0-Rik - 

Zk 

+  K4  ln(l  —  Zk)- 


K2 Zk  +  K3  In  (zk) 


The  unknown  quantities  in  the  combined  model  may  be  es¬ 
timated  using  a  system  identification  procedure.  This  model 
has  the  advantage  of  being  “linear  in  the  parameters”;  that  is, 
the  unknowns  occur  linearly  in  the  output  equation.  Given 
a  set  of  N  cell  input-output  three-tuples  {yk,  ik ,  Zk},  the  pa¬ 
rameters  may  be  solved  for  in  closed  form  using  a  result 
from  least-squares  estimation.  This  simple  off-line  (batch) 
method  is  as  follows:  We  first  form  the  vector 

Y  =  [yu  yi,  ■■■  ,  jn]T, 
and  the  matrix 


Eq.  (4)  is  the  basis  for  including  SOC  in  the  state  vector  of 
the  cell  model  as  it  is  in  state  equation  format  already,  with 
SOC  as  the  state  and  4  as  the  input. 

The  cell  model  may  be  completed  by  adding  additional 
states,  as  necessary,  and  an  output  equation.  Here,  we  first 
revisit  an  output  equation  from  an  earlier  paper  [5],  and  show 
how  it  may  be  enhanced.  Next,  we  add  a  state  to  the  model 
to  account  for  cell  hysteresis.  Thirdly,  we  add  dynamics  to 
the  state  to  model  cell  terminal  voltage  relaxation.  Finally, 
we  discuss  adding  temperature  dependence  to  the  models. 

3.3.  Models  with  only  SOC  as  a  state 

The  first  three  model  structures  that  we  investigate  have 
state  vector  Xk  =  Zk •  That  is,  the  only  state  in  the  state 
equation  (1)  is  SOC.  These  models  can  estimate  cell  terminal 
voltage  in  a  limited  way,  and  are  improved  upon  later  using 
multiple  states. 

3.3.1.  The  combined  model 

With  SOC  available  as  part  of  the  model  state,  terminal 
voltage  may  be  predicted  in  a  number  of  different  ways. 
Several  different  forms  are  adapted  from  reference  [19]. 

•  Shepherd  model:  y^  =  Eq  —  Rik  —  Ki/zk- 

•  Unnewehr  universal  model:  yk  =  E$  —  Ri ^  —  KiZk- 

•  Nernst model:  yk  =  Eo-Rijc-\-K2  ln(zk)  +  K3  ln(l  —  Zk)- 


H  =  \h\,  h2, . . .  ,  /*n]T • 

The  rows  of  H  are 

Yf,  ij ,  —,Zj,  In  (zj),  ln(l  -  z/)  , 

Zj  J 

where  ft  is  equal  to  ij  if  ij  >  0,  ij  is  equal  to  ij  if  ij  < 
0,  else  ft  and  ij  are  zero.  Then,  we  see  that  Y  =  HO, 

where  0T  =  [Ko,  R+,  R~ ,  K\,  Kl,  K3,  K4]  is  the  vector 
of  unknown  parameters.  Using  a  result  from  least-squares 
estimation  theory,  we  solve  for  the  parameters  0  using  the 
known  matrices  Y  and  H  as  0  =  (HT H)~l HTY. 

Results  comparing  the  combined  model  cell  voltage  es¬ 
timation  with  the  cell’s  true  voltage  for  the  pulsed-current 
test  are  shown  in  Fig.  6(a)  and  (d).  Fig.  6(a),  shows  the  com¬ 
parison  over  discharge  pulses,  and  Fig.  6(d)  shows  the  com¬ 
parison  over  charge  pulses.  The  general  shape  of  the  cell 
response  and  the  model  output  is  the  same,  although  many 
details  of  the  true  cell  response  are  missing.  These  will  be 
improved  with  other  models.  The  root-mean- squared  model 
estimation  error  over  the  test  shown  in  Fig.  6  is  listed  in 
Table  1 .  Results  comparing  the  combined  model  cell  voltage 
estimation  with  the  cell’s  true  voltage  for  one  cycle  of  the 
UDDS  test  is  shown  in  Fig.  8(a).  Similar  comments  apply. 

3.3.2.  The  simple  model 

With  parameter  values  fit  to  the  “combined  model”,  we 
can  evaluate  its  component  terms  for  further  insight.  The 
model  output  equation  may  be  divided  into  two  additive 
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Table  1 

Performance  of  the  model  structures  when  predicting  cell  voltage 


Model  structure 

Root-mean- squared  (RMS)  modeling  error 

Multi-rate  pulse  test  (mV) 

Multi-cycle  UDDS  test  (mV) 

— 30°C 

0°C 

25  °C 

Combined  model 

34.7 

50.1 

24.1 

23.3 

Simple  model 

36.2 

165.8 

26.9 

22.4 

Zero-state  hysteresis  model 

21.5 

62.2 

24.6 

22.3 

One-state  hysteresis  model 

21.5 

48.7 

14.1 

14.0 

Enhanced  self-correcting,  rif  =  2 

13.8 

39.0 

14.5 

7.2 

Enhanced  self-correcting,  rif  =  4 

6.7 

35.3 

4.3 

4.2 

The  root-mean-squared  value  of  modeling  error,  calculated  as  true  cell  voltage  minus  model  voltage,  and  judged  over  the  SOC  range  5-95%,  is  tabulated. 


Comparing  "combined"  OCV  with  real  OCV 


Fig.  3.  OCV  as  a  function  of  SOC,  compared  with  terms  from  the 
“combined  model.” 

parts:  one  part  depending  only  on  SOC,  and  another  depend¬ 
ing  only  on  4: 

yk  =  K0  -  —  -  K2Zk  +  K2  ln(zt)  +  K4  ln(l  -  Zk) 

Zk / 

&  (zk) 

—  Rik^  • 

fn(4) 

The  part  depending  only  on  SOC  bears  closer  examination. 
When  values  are  fit  to  parameters  {A'o, . . .  ,  K 4},  we  plot 
the  part  denoted  “fn(z^)”  versus  Zk  (Fig-  3).  Overlaid  is  the 
open-circuit  voltage  curve  as  a  function  of  SOC.2  We  see 


2  The  open-circuit  voltage  as  a  function  of  state-of-charge  for  these  cells 
as  plotted  in  Fig.  3  is  an  empirical  relationship  found  by  cell  testing.  First, 
the  cell  was  fully  charged  (constant  current  to  4.2  V,  constant  voltage 
to  200  mA).  Then,  the  cell  was  discharged  at  the  C/25  rate  until  fully 
discharged  (3.0  V).  The  cell  was  then  charged  at  the  C/25  rate  until  the 
voltage  was  4.2  V.  The  low  rates  were  used  to  minimize  the  dynamics 
excited  in  the  cells.  The  cell  voltage  as  a  function  of  state  of  charge  under 
discharge  and  under  charge  were  averaged  to  compute  the  OCV.  This  has 
the  effect  of  eliminating  to  the  greatest  extent  possible  the  presence  of 
hysteresis  and  ohmic  resistance  in  the  final  function.  For  the  purpose  of 
computations  involving  OCV,  the  final  curve  was  digitized  at  200  points 
and  stored  in  a  table.  Linear  interpolation  is  used  to  look  up  values  in 
the  table. 


Fig.  4.  Equivalent  circuit  implemented  by  “simple”  model,  and  approxi¬ 
mated  by  “combined”  model.  The  diodes  are  ideal.  R+  is  the  discharge 
resistance,  and  R~  is  the  charge  resistance. 

that  the  part  of  yk  that  is  a  function  of  SOC  is  attempting 
to  fit  the  OCV(zk)  curve.  So,  an  easier  and  more  accurate 
implementation  of  the  combined  model  is 


yk  =0  CV(zk)-Rik-  (5) 

This  output  equation  is  drawn  as  an  equivalent  circuit  in 
Fig.  4,  where  different  charge/discharge  resistances  may  be 
used.  We  call  the  model  comprised  of  Eqs.  (4)  and  (5)  the 
“simple  model”. 

This  model  type  is  also  linear  in  the  parameters.  Off-line 
system  identification  is  done  as  follows:  We  first  form  the 
vector 

Y=  [yi  —  OCV(zi),  y2  -  OCV(z2),  OCV(zn)]T, 

and  the  matrix 
H  =  [h\,  h2,  -  -  -  ,  %]T 
The  rows  of  H  are 


Again,  we  see  that  Y  =  HO,  where  0T  =  [R+,  R~]  is  the 
vector  of  unknown  parameters.  We  solve  for  the  parameters 
0  using  the  known  matrices  Y  and  H  as  0  =  (HT H)~l  HTY. 

Results  comparing  the  simple  model  cell  voltage  estima¬ 
tion  with  the  cell’s  true  voltage  for  the  pulsed-current  test 
are  shown  in  Fig.  6(b)  and  (e).  Fig.  6(b),  shows  the  com¬ 
parison  over  discharge  pulses,  and  Fig.  6(e)  shows  the  com¬ 
parison  over  charge  pulses.  The  RMS  cell  model  estimation 
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Charge  and  discharge  curves  Hysteresis  versus  SOC 


Fig.  5.  Hysteresis  curves.  Plots  of:  (a)  the  discharge/charge  curves;  (b)  the  hysteresis  level. 


error  over  the  test  shown  in  Fig.  6  is  listed  in  Table  1.  We 
see  that  the  simple  model  slightly  under-performs  the  com¬ 
bined  model  in  most  cases,  most  likely  because  the  com¬ 
bined  model  over-fits  the  data.  We  prefer  the  simple  model 
to  the  combined  model  because  we  feel  that  it  generalizes 
better,  and  because  it  is  less  complex  to  implement.  Results 
comparing  the  simple  model  cell  voltage  estimation  with  the 
cell’s  true  voltage  for  one  cycle  of  the  UDDS  test  are  shown 
in  Fig.  8(b).  Similar  comments  apply. 

3.3.3.  The  zero -state  hysteresis  model 

Examination  of  the  results  in  Fig.  6(a)  and  (b)  and  (d) 
and  (e)  exposes  some  flaws  in  the  models.  One  subtle  effect, 
which  has  serious  consequences  when  predicting  SOC,  may 
be  seen  during  rest  periods.  Following  a  discharge,  the  cell 
voltage  always  relaxes  to  a  value  less  than  the  true  OCV 
for  that  SOC,  and  following  a  charge,  the  cell  voltage  al¬ 
ways  relaxes  to  a  value  greater  than  the  true  OCV.  This  is 
explained  by  a  hysteresis  effect  occurring  in  the  cell  that  is 
not  modeled  in  the  combined  or  simple  models.3 

The  term  hysteresis  is  derived  from  the  Greek  hustereia 
“to  arrive  late”.  The  cell  voltage  lags  the  predicted  voltage 
in  some  sense.  It  may  also  be  defined  as  a  characteristic  of  a 
system  in  which  a  change  in  the  direction  of  the  independent 
variable  leads  to  the  dependent  variable  failing  to  retrace  the 
path  it  passed  in  the  forward  direction.  (For  a  good  paper 
describing  electro-chemical  hysteresis,  see  Ref.  [20].) 

We  illustrate  this  effect  by  showing  the  charge/discharge 
curves  at  the  C/25  rate  and  room  temperature  in  Fig.  5(a). 
The  terminal  voltage  for  discharge  is  the  lower  curve;  for 
charge  the  upper  curve.  Two  different  terminal  voltages  exist 
at  each  SOC.  Half  the  difference  between  these  voltages  is 


3  We  note  in  passing  that  hysteresis  is  not  a  phenomenon  generally 
associated  with  lithium-ion  systems,  since  most  applications  have  been  in 
the  light  portable  electronics  area  where  SOC  accuracy  is  not  as  critical 
as  in  the  HEV  application  and  where  temperatures  are  not  as  extreme.  It 
is,  however,  very  pronounced  at  low  temperatures  and  can  lead  to  SOC 
errors  as  large  as  ±40%  if  the  estimate  is  based  simply  on  OCV  (even 
with  full  cell  relaxation).  The  reason  is  the  spread  between  the  charge 
and  discharge  characteristics  coupled  with  the  flat  nature  of  the  curves 
between  10  and  90%  SOC. 


the  polarization  voltage  of  the  cell.  Only  a  small  part  of  the 
polarization  is  due  to  Rik  drop  (about  2.5  mV  here)  and  the 
remainder  is  due  to  hysteresis  effects.  The  hysteresis  level, 
with  Rik  subtracted  out,  is  plotted  versus  SOC  in  Fig.  5(b). 
We  have  found  that  cell  voltage  hysteresis  is  considerably 
larger  at  low  temperatures,  and  that  we  must  include  hys¬ 
teresis  in  our  cell  model  for  good  SOC  estimation. 

These  curves  comprise  the  major  hysteresis  loop ,  corre¬ 
sponding  to  full  cell  charge  and  discharge.  Minor  hysteresis 
loops  are  encountered  when  a  partial  charge  is  followed  by  a 
partial  discharge,  and  vice  versa.  The  polarization  does  not 
immediately  flip  sign  upon  a  current  reversal,  but  slowly  de¬ 
cays  from  one  leg  of  the  major  hysteresis  loop  to  the  other.4 
To  capture  the  dynamics  of  the  gradual  decay  in  voltage  we 
need  to  add  one  or  more  states  to  the  model,  which  is  done 
in  Section  3.4.1.  Here,  we  explore  a  simpler  version  as  a 
proof-of-concept.  This  model  adds  no  additional  states  for 
hysteresis,  so  is  named  the  “zero-state  hysteresis  model”. 

A  basic  model  of  hysteresis  simply  adds  a  term  to  the 
output  equation 


yk  =  OCV(z*)  -  sjcM(zk)  -  Rik , 

where  Sk  represents  the  sign  of  the  current  (with  memory 
during  a  rest  period).  For  some  s  sufficiently  small  and  pos- 


itive, 

1, 

ik  >  £, 

Sk  =  ‘ 

-1, 

77' 

A 

1 

00 

Sk- 1, 

141  <  £• 

M(zk)  is  half  the  difference  between  the  two  legs  of  the 
charge/discharge  curve,  minus  the  Rik  loss,  and  is  plotted  in 
Fig.  5(b).  Here,  we  use  a  constant  value  for  M. 


4  Successively  smaller  concentric  minor  loops  may  be  obtained  by 
alternating  shorter  and  shorter  charge  and  discharge  pulses,  eventually 
converging  on  the  mean  of  the  two  values  of  Fig.  5(a)  at  each  SOC. 
Therefore,  we  compute  OCV  as  a  function  of  SOC  as  the  mean  of  the 
two  legs  of  the  major  hysteresis  loop. 
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This  model  type  is  also  linear  in  the  parameters.  Off-line 
system  identification  is  done  as  follows:  We  first  form  the 
vector 

Y=  [VI  — OCV(zi),  V2  -  OCV(z2),  •  •  •  ,  Vn  -  OCV(zN)]T, 
and  the  matrix 

H  =  [huh!,...  ,/iN]T. 

The  rows  of  H  are 


Modeling  discharge:  Combined  model 


Modeling  discharge:  Simple  model 


Modeling  discharge:  Zero  state  hysteresis  model 


Again,  we  see  that  Y  =  HO,  where  0T  =  [7?+,  R~ ,  M]  is  the 
vector  of  unknown  parameters.  We  solve  for  the  parameters 
0  using  the  known  matrices  Y  and  H  as  0  =  (HT H)~l  HTY. 

Results  comparing  the  zero-state  hysteresis  model  cell 
voltage  estimation  with  the  cell’s  true  voltage  for  the 
pulsed-current  test  are  shown  in  Figs.  6(c)  and  (f).  Fig.  6(c), 
shows  the  comparison  over  discharge  pulses,  and  Fig.  6(f) 
shows  the  comparison  over  charge  pulses.  The  RMS  cell 
model  estimation  error  over  the  test  shown  in  Fig.  6  is 
listed  in  Table  1.  Performance  of  the  zero-state  hysteresis 
model  is  consistently  better  than  that  of  the  simple  model. 
Results  comparing  the  zero-state  hysteresis  model  cell  volt¬ 
age  estimation  with  the  cell’s  true  voltage  for  one  cycle  of 


Modeling  charge:  Combined  model 


4.2 

4.0 

—  3.8 
> 

ji=i 

CD 

I53'6 

o 

>  3.4 

3.2 

3.0 

True  voltage 
-  Estimated  voltage  : 

0  10  20  30  40  50  60  70  80  90 

(d)  Time  (min.) 


Modeling  charge:  Simple  model 


Modeling  charge:  Zero  state  hysteresis  model 


Fig.  6.  Results  of  cell  modeling  using  models  with  only  SOC  as  a  state  for  the  pulsed-current  cell  tests.  Discharge  portion  of  test  is  shown  in  (a)-(c); 
charge  portion  of  test  is  shown  in  (d)-(f).  The  gray  line  is  the  measured  cell  voltage,  and  the  black  line  is  the  model  prediction. 
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the  UDDS  test  are  shown  in  Fig.  8(c).  Similar  comments 
apply. 

3.4.  Models  with  SOC  and  additional  states 


In  order  to  better  estimate  cell  voltage  effects  that  are  cou¬ 
pled  to  the  history  of  the  cell’s  input  current,  we  must  make 
modifications  to  the  model  state  equation  (1).  We  examine 
two  additions  in  the  following  sections. 


3.4.1.  The  one-state  hysteresis  model 

The  zero-state  hysteresis  model  is  an  improvement  over 
the  simple  model,  but  only  crudely  approximates  the  under¬ 
lying  phenomenon.  Whereas  the  level  of  hysteresis  slowly 
changes  as  the  cell  is  charged  or  discharged,  the  model  es¬ 
timates  hysteresis  as  immediately  flipping  between  its  max¬ 
imum  positive  and  negative  values  when  the  sign  of  current 
changes. 

The  slow  transition  may  be  modeled  by  adding  a  “hys¬ 
teresis  state”  to  the  model  state  equation  (1).  The  hystere¬ 
sis  state  is  not  a  differential  equation  in  time,  but  in  SOC 
(or,  ampere-hours).  Let  h(z,t)  be  the  hysteresis  voltage  as  a 
function  of  SOC  and  time,  and  let  z  =  dz/dt.  Then, 

d h(z  t) 

— —3—  =  y  sgn(z)(M(z,  z)  -  h(z ,  t))9 

d  z 

where  M(z,  z)  is  a  function  that  gives  the  maximum  po¬ 
larization  due  to  hysteresis  as  a  function  of  SOC  and  the 
rate-of-change  of  SOC.  Specifically,  M(z,  z)  is  positive  for 
charge  (z  >  0)  and  is  negative  for  discharge  (z  <  0).  The 
M(z,  z)  —h(z,  t)  term  in  the  differential  equation  states  that 
the  rate-of-change  of  hysteresis  voltage  is  proportional  to  the 
distance  away  from  the  major  hysteresis  loop,  leading  to  a 
kind  of  exponential  decay  of  voltage  to  the  major  loop.  The 
term  in  front  of  this  has  a  positive  constant  y,  which  tunes 
the  rate  of  decay,  and  sgn(z),  which  forces  the  equation  to 
be  stable  for  both  charge  and  discharge. 

In  order  to  fit  the  differential  equation  for  h(z,t)  into  our 
model,  we  must  manipulate  it  to  be  a  differential  equation  in 
time,  not  in  SOC.  We  accomplish  this  by  multiplying  both 
sides  of  the  equation  by  dz/dt 


d h(z,  t)  d z 
d  z  d  t 


ysgn(z)(M(z,z)  -  h(z,t )) 


d  z 
dt 


Note  that  dz/dt  =  —  rjii(t)/Cn,  and  that  z  sgn(z)  =  |z|.Thus, 


Note  that  this  is  a  linear-time-varying  system  as  the  factors 
multiplying  the  state  and  input  change  with  ik  and  hence 
with  time.  If  we  define  F(ik)  =  exp(—  \rjiikyAt/Cn\),  then 
the  overall  state-space  equations  for  the  one-state  hysteresis 
model  are 


h 

M(z ,  z)  J  ’ 

yk  =  OCVfeO  -  Rik  +  hk. 

Results  comparing  the  one- state  hysteresis  model  cell 
voltage  estimation  with  the  cell’s  true  voltage  for  the 
pulsed-current  test  are  shown  in  Fig.  7(a)  and  (d).  Fig.  7(a), 
shows  the  comparison  over  discharge  pulses,  and  Fig.  7(d) 
shows  the  comparison  over  charge  pulses.  The  RMS  cell 
model  estimation  error  over  the  test  shown  in  Fig.  7  is 
listed  in  Table  1.  Performance  of  the  one-state  hysteresis 
model  is  consistently  better  than  the  simpler  models.  Re¬ 
sults  comparing  the  one-state  hysteresis  model  cell  voltage 
estimation  with  the  cell’s  true  voltage  for  one  cycle  of  the 
UDDS  test  are  shown  in  Fig.  8(d).  Similar  comments  apply. 

3.4.2.  The  enhanced  self -correcting  (ESC)  model 

A  significant  element  missing  from  these  models  is  a  de¬ 
scription  of  time  constants  during  pulsed  current  events.  If 
a  cell  is  allowed  to  rest,  it  takes  some  time  for  the  voltage 
to  completely  relax  to  its  rest  voltage.  If  a  cell  is  pulsed 
with  current,  it  takes  time  for  the  voltage  to  converge  to  its 
steady-state  level.  These  time  constants,  which  describe  the 
phenomenon  we  henceforth  refer  to  as  the  relaxation  effect, 
may  be  implemented  as  a  low-pass  filter  on  ik.  Since  the 
cell  model  must  accurately  predict  its  behavior  in  a  dynamic 
HEV  environment,  we  find  it  is  essential  to  include  relax¬ 
ation  effects. 

Early  attempts  to  model  the  relaxation  effect  included 
filtering  the  state-of-charge  as  well  as  the  input  current  (cf. 
the  “filter  state”  cell  model  [5]).  While  this  model  could 
fit  voltage  data  reasonably  well,  it  had  the  unfortunate  side 
effect  that  SOC  estimation  using  an  EKF  was  not  reliable. 
The  output  equation  had  the  form 


hk+ 1 
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yk  =  OC Y(zk)  +  Alt (zk)  +  filt(/*)  -  Rk , 


This  may  be  converted  into  a  difference  equation  for  our 
discrete-time  application  using  standard  techniques  (assum¬ 
ing  that  i(t)  and  M(z,  z)  are  constant  over  the  sample  period): 


hk+ 1 


^  hk 
mhyAt 


M(z,  z). 


where  filt(-)  is  some  dynamic  operation  filtering  its  operand. 
The  EKF  had  problems  with  this  model  because  the 
credit/blame  portion  of  the  algorithm  could  not  reliably  de¬ 
termine  whether  model  error  was  due  to  bad  SOC  through 
the  OCV(-)  function  or  through  the  filt(z&)  function,  or 
due  to  bad  filter  states  in  the  filt(4)  function.  In  particular, 
some  simple  cell  tests  running  the  EKF  showed  a  lack  of 
robustness: 
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Modeling  discharge:  One  state  hysteresis  Modeling  charge:  One  state  hysteresis 


0  10  20  30  40  50  60  70  0  10  20  30  40  50  60  70  80  90 

(a)  Time  (min.)  (d)  Time  (min.) 


Modeling  discharge:  ESC,  2  filter  states  Modeling  charge:  ESC,  2  filter  states 


0  10  20  30  40  50  60  70  0  10  20  30  40  50  60  70  80  90 

(b)  Time  (min.)  (e)  Time  (min.) 


Modeling  discharge:  ESC,  4  filter  states  Modeling  charge:  ESC,  4  filter  states 
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(c)  Time  (min.)  (f)  Time  (min.) 

Fig.  7.  Results  of  cell  modeling  using  models  with  multiple  states  for  the  pulsed-current  cell  tests.  Discharge  portion  of  test  is  shown  in  (a)-(c);  charge 
portion  of  test  is  shown  in  (d)-(f).  The  gray  line  is  the  measured  cell  voltage,  and  the  black  line  is  the  model  prediction. 


1.  A  constant-current  discharge/charge  should  make  the 
SOC  ramp  down/up  at  the  slope  i /  Cn  (A/Ah).  In  practice, 
the  slope  using  the  filter  state  model  was  often  wrong. 

2.  During  a  rest  period,  cell  terminal  voltage  converges  to 
OCV  (neglecting  hysteresis  effects)  and  estimated  SOC 
should  converge  to  the  SOC  predicted  by  OCV.  In  the 
implementation,  we  observed  SOC  to  drift  considerably, 
not  converging  to  the  correct  value. 

The  model  that  we  will  develop  in  this  section,  called  the 
“enhanced  self-correcting  model”,  forces  y&  to  converge  to 
OCV  after  a  rest  period  and  it  forces  y \  to  converge  to 


OCV  —  Rik  for  a  constant-current  discharge/charge.  To  meet 
these  requirements,  with  hysteresis  added,  the  output  equa¬ 
tion  needs  to  have  the  form 

yk  =  OCV (zk)  +  Jik^  +  filt(4)  -  Rik  ■ 

fn  (Zk)  fn  (zkJk)  fn  Ok) 

In  this  equation,  SOC  and  hysteresis  contribute  the  long-term 
dc  level  (bias)  to  the  output  and  4  and  its  history  contribute 
the  short-term  variation  around  this  level.  SOC  itself  is  no 
longer  filtered  as  in  the  “filter  state  model” — it  makes  no 
sense  to  have  a  moving  bias  point. 
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Modeling  (zoom):  Combined  model  Modeling  (zoom):  Simple  model 


Modeling  (zoom):  Zero  state  hysteresis  model  Modeling  (zoom):  One  state  hysteresis 


Modeling  (zoom):  ESC,  2  filter  states  Modeling  (zoom):  ESC,  4  filter  states 
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Fig.  8.  Results  of  cell  modeling  for  the  UDDS  cell  tests.  The  gray  line  is  the  measured  cell  voltage,  and  the  black  line  is  the  model  prediction. 


The  filter  “filt(-)”  must  satisfy  two  criteria:  (1)  after  a 
long  rest  period  its  output  must  be  zero  so  that  y \  — > 
OCV  +  hk\  (2)  during  a  constant-current  discharge/charge, 
its  output  must  converge  to  zero  so  that  yk  — >  OCV  + 
hk  —  Rik •  The  first  criterion  is  satisfied  by  a  stable  linear 
filter,  and  the  second  is  satisfied  by  a  linear  filter  with  zero 
dc  gain.  Both  of  these  may  be  enforced  in  the  filter  de¬ 
sign. 

A  linear  filter  may  be  implemented  in  a  state- space  form 
as 

fk+ 1  =  Mfk  +  Bfik, 

y{  =  G/*> 


where  Af  is  the  state-transition  matrix  of  the  filter,  Bf  is  the 
input  matrix  of  the  filter,  G  is  the  output  matrix  of  the  filter, 
and  fk  is  the  filter  state.  The  eigenvalues  of  the  Af  matrix 
are  the  “poles”  of  the  filter  and  determine  its  stability.  The 
filter  is  stable  if  max|eig(Af)  |  <  1.  The  location  of  poles  de¬ 
termine  the  system’ s  dynamic  behavior:  poles  near  + 1  have 
slowly  decaying  dynamics,  poles  near  zero  decay  quickly, 
negative  poles  oscillate.  Complex-conjugate  poles  also  oscil¬ 
late,  and  do  not  appear  to  improve  performance  by  their  in¬ 
clusion.  Therefore,  it  is  sufficient  to  have  an  Af  matrix  of  the 
form  Af  =  diag((z),  where  a  is  a  vector  comprising  the  pole 
locations.  Stability  is  ensured  if  all  —1  <  glj  <  1.  The  Bf 
matrix  may  be  chosen  arbitrarily  so  long  as  no  entry  is  zero. 
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Now  that  we  have  guaranteed  stability,  y&  OCV  +  hk 
during  rest.  By  carefully  selecting  the  G  matrix,  we  can 
ensure  a  zero  dc-gain  of  the  filter  so  that  y&  OCV  +  hk  — 
Rik  during  constant-current  profiles.  The  gain  of  the  filter  is 
(where  nf  is  the  number  of  filter  states,  determined  during 
system  identification  to  fit  cell  data): 


states.  We  do  not  see  much  improvement  by  increasing  nf  be¬ 
yond  4.  Results  comparing  the  ESC  model  cell  voltage  esti¬ 
mation  with  the  cell’ s  true  voltage  for  one  cycle  of  the  UDDS 
test  are  shown  in  Figs.  8(e)  for  rif  =  2  and  (f)  for  rif  =  4. 

3.5.  Adding  temperature  dependence  to  the  models 
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If  we  let  g\  through  gnf-i  be  found  by  a  system-identification 
procedure  and  assuming  that  Bf  =  [  1  •  •  •  1  ]T,  then  the 

zero  dc-gain  constraint  fixes  gnf  as 
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So,  the  full  self-correcting  model  is 
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Results  comparing  the  ESC  model  cell  voltage  estimation 
with  the  cell’s  true  voltage  for  the  pulsed-current  test  are 
shown  in  Figs.  7(b)  and  (e)  for  nf  =  2  and  in  Figs.  7(c) 
and  (f)  for  nf  =  4.  Figs.  7(b)  and  (c),  show  the  comparison 
over  discharge  pulses,  and  Figs.  7(e)  and  (f)  show  the  com¬ 
parison  over  charge  pulses.  The  RMS  cell  model  estimation 
error  over  the  tests  shown  in  Fig.  7  are  listed  in  Table  1.  Per¬ 
formance  is  significantly  improved  by  the  addition  of  filter 


Thus  far,  we  have  discussed  only  how  to  model  cell  dy¬ 
namics  at  one  specific  temperature.  We  now  embark  on  a 
brief  discussion  on  how  to  incorporate  temperature  depen¬ 
dence  into  the  models. 

A  very  simple  method,  and  the  one  we  tried  first,  was  to 
use  a  table  of  different  models,  where  each  model  had  pa¬ 
rameters  optimized  for  a  specific  temperature.  For  example, 
we  used  sixteen  models  over  the  temperature  range  —30  to 
45  °C  in  increments  of  5  °C.  This  worked  well  so  long  as 
the  cell  under  test  had  temperature  equivalent  to  one  of  the 
sixteen  stored  models.  If  the  temperature  was  between  two 
stored  model  values,  we  linearly  interpolated  model  param¬ 
eters  between  the  parameters  of  the  models  in  the  table.  This 
did  not  work  well. 

We  found  that  two  adjacent  models  in  the  table  did  not 
necessarily  have  similar  parameters.  Individually  optimiz¬ 
ing  model  parameters  at  specific  temperatures  resulted  in 
values  that  were  over-fit  to  the  data  and  did  not  generalize 
well  to  cases  not  previously  seen.  We  remedied  the  problem 
by  performing  joint  optimization  over  the  entire  temperature 
range,  where  every  parameter  was  represented  by  a  contin¬ 
uous  polynomial  of  temperature  (fourth  order).  This  forced 
nearby  models  to  have  similar  parameter  values.  Although 
joint  optimization  did  not  result  in  modeling  errors  as  low 
as  when  individually  optimized,  the  generalization  perfor¬ 
mance  was  much  better. 

Fig.  9  shows  the  results  for  joint  optimization  over 
the  entire  temperature  range,  and  Table  1  lists  some  nu¬ 
meric  values  corresponding  to  the  plots.  The  cell  data  was 
collected  from  UDDS  tests  similar  to  those  described  in 
Section  3.1,  but  performed  at  16  controlled  temperatures 
from  —30  to  45  °C,  in  increments  of  5° C.  At  lower  temper¬ 
atures,  the  magnitude  of  the  current  had  to  be  scaled  down 
so  as  not  to  exceed  voltage  limits  (due  to  increased  cell 


RMS  modeling  error  versus  temperature 


Summary  modeling  performance  over  temperature 


(b)  Average  RMS  modeling  error  (mV) 


Fig.  9.  Results  when  modeling  over  a  temperature  range.  Frame  (a)  shows  the  individual  modeling  results,  and  frame  (b)  compares  the  average  modeling 
results. 
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resistance  at  lower  temperatures),  and  hence  more  UDDS 
cycles  had  to  be  completed  to  cover  the  desired  SOC  range, 
but  the  tests  were  otherwise  the  same.  In  frame  (a),  the 
RMS  modeling  errors  for  the  jointly  optimized  models  are 
plotted  versus  temperature  for  the  different  cell  models. 
In  frame  (b),  the  average  RMS  modeling  errors  over  all 
temperatures  are  presented  as  a  bar  graph. 

We  see  that  the  combined  model  appears  to  perform  well. 
However,  this  is  an  artifact.  The  {Ko,  •  •  •  ,  K4}  points  were 
found  to  over-fit  the  measured  data,  so  that  the  resulting 
curve  plotted  with  their  values  did  not  resemble  OCV  to  any 
degree  of  fidelity.  While  the  simple  model  has  worse  nu¬ 
meric  indicators  of  performance,  it  generalized  better  than 
the  combined  model.  By  adding  the  crude  model  of  hystere¬ 
sis,  we  see  a  significant  performance  jump,  especially  at  cold 
temperatures  where  hysteresis  is  most  evident  in  our  cells’ 
dynamics.  Adding  the  dynamics  of  a  state  to  the  hysteresis 
model  also  improves  performance,  again  with  the  greatest 
gains  at  low  temperatures.  Filter  states  also  contribute  to  per¬ 
formance  gains.  The  final  model,  enhanced-self-correcting 
with  rif  =  4  gave  the  best  performance  in  all  cases.  By 
increasing  the  number  of  filter  states  we  would  expect  con¬ 
tinued  performance  gains,  at  the  cost  of  greater  complexity. 
Note  that  the  cold-temperature  performance  is  not  improved 
as  much,  in  relative  terms,  by  adding  filter  states,  so  it  is 
likely  that  some  cold-temperature  phenomena  is  not  yet 
being  modeled  well.  Conceivably,  a  second  hysteresis  state 
could  be  added  to  the  model  to  improve  performance  here. 
We  have  not  investigated  this  possibility  as  yet. 


4.  System  identification 

The  first  three  system  models  introduced  in  this  paper  are 
“linear  in  the  parameters”.  This  makes  identifying  the  values 
of  the  model  parameters  straightforward  using  least-squares 
estimation,  and  has  been  discussed  earlier.  When  the  model 
is  not  linear  in  the  parameters,  as  in  the  remaining  system 
models,  this  method  may  not  be  used.  We  must  turn  to  more 
advanced  methods. 

Here  we  look  at  one  method  in  particular.  We  know  that 
a  Kalman  filter  or  extended  Kalman  filter  may  be  used  to 
estimate  the  state  of  a  dynamic  system  given  noisy  mea¬ 
surements;  e.g.,  to  estimate  the  cell  SOC.  We  may  also  use 
an  extended  Kalman  filter  to  perform  system  identification 
given  clean  measurements.  To  do  so,  we  require  a  state- space 
model  describing  the  dynamics  of  the  parameters  0  of  the 
system  model.  We  will  use  the  Kalman  filter  as  an  optimum 
observer  of  these  parameter  values,  creating  an  estimate  0. 
In  electro-chemical  cells,  the  true  parameters  will  change 
only  very  slowly,  so  we  model  them  as  constant  with  some 
small  perturbation: 

0k+ 1  =0k  +  rk- 

The  small  white  noise  input  r k  is  fictitious,  but  models 
the  slow  change  in  the  parameters  of  the  system  plus  the 


Table  2 

Summary  of  the  nonlinear  extended  Kalman  filter  for  system  identification 
[21] 

Nonlinear  state- space  modela 


Ok+i  =0k  +  rk 

dk  =  g(xk,  uk,Ok)  +  ek 


Definition 

_  d g(xk,uk,0) 
Lk  ~  dO 

Initialization 
For  k  =  0,  set 


et  =  E[0O] 

=  E[(0O  -  eo+)(0o  -  e+)T] 
Computation 

For  k  =  1,2,...  compute 

State  estimate  time  update:  Of  =  0k_x 


Error  covariance  time  update:  E~  =  Ef  +  £r 

6,k  6,k—  1 

Kalman  gain  matrix:  L\  =  rrjt(Cf)T[CfX'--i(Cf)T  +  Se]~ 1 

State  estimate  measurement  update: 

0*  =  Of  +  Lk[yk  -  g(xk,  uk,  Of)] 

Error  covariance  measurement  update:  27^  =  (/  —  LekCek)Ef  ^ 

a  rk  and  ek  are  independent,  zero-mean,  Gaussian  noise  processes  of 
covariance  matrices  Er  and  Ee,  respectively. 


infidelity  of  the  model  structure  to  capture  all  of  the  cell 
dynamics. 

The  output  equation  required  for  Kalman-filter  system 
identification  must  be  a  measurable  function  of  the  system 
parameters.  We  use 

dk  =  8(%ki  Uk,  Ok)  ~\~  6ki 


where  g(-,  •,  •)  is  the  output  equation  of  the  system  model 
being  identified,  and  ek  models  the  sensor  noise  and  model¬ 
ing  error.  We  compare  <4  computed  using  0k  to  the  measured 
cell  output,  and  adapt  Ok  to  minimize  the  difference. 

We  can  create  an  extended  Kalman  filter  using  this 
state- space  model  and  cell  data  to  estimate  the  system  pa¬ 
rameters  as  summarized  in  Table  2.  We  initialize  the  state 
estimate  with  our  best  information  re.  the  state  value:  0q  = 
E[0oL  and  the  state  estimation  error  covariance  matrix: 
sf  =  E[(0  -  e+){6  -  3+)T]. 

The  time  update  propagates  the  state  estimate  as  0^  = 
0k-\  since  the  parameters  are  assumed  constant,  and  the 

error  covariance  as  E~1  =  +  Er  to  account  for  the 

0,k  0,k—  1 

added  uncertainty  due  to  the  fictitious  noise  input  The 
effect  of  adding  Ur  is  to  increase  the  estimate’s  uncertainty, 
and  to  allow  adaptation  to  0. 

The  extended  Kalman  filter  gain  matrix  is  computed  by 
linearizing  the  state-space  model’s  output  equation.  We  com¬ 
pute 


q  = 

T0 


dg(xk,  uk,0) 


d  e 


0=0,7 


ET(Ci)l[ClST(CW+Ee]- 
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Next,  we  measure  the  true  cell  output  yk  and  compare  it  to  the 
model  output  g(xk,  uk,0f).  To  do  so,  we  need  to  simulate 
the  model  in  parallel  with  the  real  cell  to  have  an  appropriate 
value  of  Xk  available.  The  difference  between  cell  output 
and  model  output  can  be  attributed  to  noises  and  modeling 
error.  The  extended  Kalman  filter  adapts  0  to  minimize  this 
difference 

ft  —  ft  T  Lkiyk  g{xk,  Uk’  ft  )]• 

Finally,  the  measurement  update  is  applied  to  the  state  error 
covariance  matrix 


o,k 


=  (/  -  nci)z: 


v— 

J0,k' 


This  has  the  effect  of  decreasing  the  modeling  uncertainty 
due  to  the  measurement  update.  All  terms  are  accounted  for, 
and  the  algorithm  is  complete. 


4.1.  Extended  Kalman-filter  system  identification  for  the 
one-state  hysteresis  model 


The  details  for  applying  this  method  to  any  particular  cell 
model  are  differentiated  by  the  calculation  of  C6k.  For  the 
one-state  hysteresis  model,  let  the  vector  of  parameters  be 

0  =  [R+  R~  M  y  ]T. 

R+  is  the  cell  resistance  when  current  is  positive,  and  R~ 
is  the  cell  resistance  when  current  is  negative.  M  is  the 
maximum  hysteresis  voltage,  and  y  is  the  hysteresis  rate 
constant,  which  is  part  of  F(ik).  To  calculate  Cek  we  require 

d  g(xk,uk,6)  dg(xk,  uk,  0)  dg(xk,  uk,  0)  dxk  ... 

-  = - - ,  (6) 

d  0  dO  dxk  dO 

d xk  _  df(xk- 1,  uk- 1,  0)  df(xk-i,uk-i,  0)  dxk-i 
d  0  ~  dO  dxk-i  d  O'  (  j 

The  derivative  calculations  are  recursive  in  nature,  and 
evolve  over  time  as  the  state  evolves.  The  term  dxo/dO  is 
initialized  to  zero  unless  side  information  gives  a  better  es¬ 
timate  of  its  value.  We  see  that  in  order  to  calculate  Ck  for 
any  specific  model  structure,  we  require  methods  to  calcu¬ 
late  the  partial  derivatives  in  (6)  and  (7).  For  the  one- state 
hysteresis  model  we  have 


4.2.  Extended  Kalman-filter  system  identification  for  the 
enhanced  self- correcting  model 


For  the  enhanced  self-correcting  model,  let  the  vector  of 
parameters  be 

6=[R+  R~  gl  •••  gnf- 1  Pi  ■  ■  ■  Pnf  M  Y  ]T, 


where  ft  =  tanh  (a)  and  a  is  the  vector  of  filter  pole  lo¬ 
cations.  We  use  the  tanh  (•)  function  during  system  identi¬ 
fication  because  it  forces  filter  poles  to  remain  within  ±1 
(i.e.,  stable)  regardless  of  the  value  of  f>.  When  calculating 
the  partial  derivatives  we  must  remember  that  since  gnf  = 

—  YTi= l  gi(l  —  °7if)/(l  ~  oti),  it  is  not  independently  iden¬ 
tified  but  is  computed  from  the  g\,  . . .  ,  gnf- 1  terms.  This 
also  forces  the  derivatives  to  be  more  complicated  than  a 
quick  glance  would  indicate.  That  is, 


dg(xk,  uk ,  0) 
dgl 


=  fk,t 


1  ~  anf 

1  —  Oil 


fk , 


rif  •> 


and  so  forth.  Also,  since  oti  =  tanh  (ft),  it  can  never  be 
unity,  so  division  by  zero  is  impossible  in  the  derivative 
computation.  With  this  in  mind,  the  partial  derivatives 
in  (6)  and  (7)  may  be  computed  in  a  straightforward 
way. 


5.  Conclusions 

This  paper  has  proposed  five  mathematical  state-space 
structures  for  the  purpose  of  modeling  LiPB  HEV  cell  dy¬ 
namics  for  their  eventual  role  in  HEV  BMS  algorithms. 
Models  with  a  single-state  are  very  simple,  but  perform  the 
poorest.  Adding  hysteresis  and  filter  states  to  the  model  aids 
performance,  at  some  cost  in  complexity. 

We  have  also  seen  how  to  identify  the  parameters  of  the 
cell  models  given  cell-test  data.  Models  that  are  “linear  in  the 
parameters”  may  have  their  parameters  fit  in  a  very  straight¬ 
forward  way  using  methods  from  least- squares-estimation 
theory.  Models  with  more  dynamics  than  simply  SOC  re¬ 
quire  more  sophisticated  techniques.  One  possibility  is  to 
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Note  that  the  dOCV(zk)/8zk  term  is  never  needed,  as  it 
always  multiplies  zero.  For  this  particular  model,  we  can 
simplify  the  calculations  by  removing  the  multiplies  by  zero: 


use  an  extended  Kalman  filter  to  identify  the  cell  parameters 
in  an  on-line  or  off-line  manner. 
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In  the  third  paper  [4],  we  will  employ  extended  Kalman 
filtering  from  [3],  using  the  cell  models  developed  here,  to 
implement  HEV  BMS  algorithms.  We  will  see  how  to  use 
EKF  to  estimate  SOC  and  all  other  model  states  as  the  sys¬ 
tem  operates.  This  model  state  will  then  allow  us  to  accu¬ 
rately  compute  a  dynamic  estimate  of  available  power.  We 
can  additionally  employ  a  technique  called  dual  extended 
Kalman  filtering  to  simultaneously  estimate  cell  state  and 
parameters,  allowing  tracking  of  cell  power  fade  and  capac¬ 
ity  fade,  for  example.  Finally,  the  parameter  data  and  SOC 
estimate  may  be  combined  to  determine  which  cells  in  the 
pack  must  have  their  charge  levels  modified  in  order  to  bring 
the  pack  into  equalization. 
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